Fly casting with ligand sliding and orientational selection supporting complex formation of a GPCR and a middle sized flexible molecule

A GA-guided multidimensional virtual-system coupled molecular dynamics (GA-mD-VcMD) simulation was conducted to elucidate binding mechanisms of a middle-sized flexible molecule, bosentan, to a GPCR protein, human endothelin receptor type B (hETB). GA-mD-VcMD is a generalized ensemble method that produces a free-energy landscape of the ligand-receptor binding by searching large-scale motions accompanied with stable maintenance of the fragile cell-membrane structure. All molecular components (bosentan, hETB, membrane, and solvent) were represented with an all-atom model. Then sampling was conducted from conformations where bosentan was distant from the binding site in the hETB binding pocket. The deepest basin in the resultant free-energy landscape was assigned to native-like complex conformation. The following binding mechanism was inferred. First, bosentan fluctuating randomly in solution is captured using a tip region of the flexible N-terminal tail of hETB via nonspecific attractive interactions (fly casting). Bosentan then slides occasionally from the tip to the root of the N-terminal tail (ligand–sliding). During this sliding, bosentan passes the gate of the binding pocket from outside to inside of the pocket with an accompanying rapid reduction of the molecular orientational variety of bosentan (orientational selection). Last, in the pocket, ligand–receptor attractive native contacts are formed. Eventually, the native-like complex is completed. The bosentan-captured conformations by the tip-region and root-region of the N-terminal tail correspond to two basins in the free-energy landscape. The ligand-sliding corresponds to overcoming of a free-energy barrier between the basins.

. Notation/terms specific to GA-mD-VcMD This table explains notation/terms that are specific GA-mD-vcMD sampling method or the current molecular system studied. Sampling method used for the current study. This is a generalized ensemble method introduced to enhance conformational sampling of a biological system, which consists of receptors and ligands in this paper.

Reaction coordinate (RC)
Given two atom groups introduced in the molecular system by a user, a reaction coordinate (RC) is defined by the inter-centroid distance between the two atom groups (see Figure S1).
In GA-mD-VcMD, conformational sampling is enhanced in a reaction-coordinate space.
Definition of an atom group is arbitrary in theory. In this study, we introduce three RCs.
Conformational sampling is enhanced in the 3D RC space. Actual definition of the three RCs is shown in Figure 1. This quantity is an ensemble of ten conformations in the 3D RC space obtained from separation simulations. These conformations are the initial conformations of diffusion simulations in the 3D RC space. Because the ten conformations are close to one another, the ten spheres look like a single conformation in Figure S6.  Figure   S4.

Zone
The 3D RC space is divided into small rectangles, and each of which is named "zone" or "virtual zone". The size of the zone may depend on its position in the 3D RC space. Actual setting of zones is given in Table S2.
Human endothelin receptor type B (hETB) A GPCR protein, which act as receptor of ligand, bosentan. This molecule has a deep pocket, to which ligands bind.

Bosentan
Ligand that binds to the deep binding pocket of hETB. This molecule is a drug molecule.

Tail
The N terminal tail of hETB, whose sequence is: Ser80, Pro81, Pro82, Arg83, Thr84, Ile85, Ser86, Pro87, Pro88, Pro89. The tail is intrinsically disordered. The N-terminal residue Ser80 is capped by an acetyl group Ace. A residue number of 79 may be assigned to Ace if necessarily (see Figure 5 for instance).

Trans-membrane helix (TM) helix region of hETB used for superposition
This region is defined in red-colored segments in Figure S7a. This region is used for conformational superposition between a snapshot and the native complex structure before computing various quantities in this study.
./012 and 3042 is a quantity to show the positional, orientational, and conformational difference of bosentan between a snapshot and the native complex. We defined two types of .
./012 is calculated using all heavy atoms of bosentan between a snapshot and the native complex structure. 3042 is done using heavy atoms in the bosentan's core region. Definition of the core region is given in Figure 3a. Note that the conformational superposition between a snapshot and the native complex structure is done for the Transmembrane helix (TM) helix region of hETB (see Figure S7a). Probability distribution function of detecting the system in a zone that involves the position ( (") , ($) , (%) ) in the 3D RC. This distribution is an equilibrated distribution at 300 K (simulation temperature). Figure S8 is the actual distribution function obtained from GA-mD-VcMD. ; Thermodynamic weight (statistical weight) assigned to snapshot at simulation temperature (300 K). This quantity is used to calculate a probability distribution function such as <=> ( 3?>2 ) or <=@ ( 3?>2 ), which are explained below. Thermal average of quantities is also calculated using ; . Cube and its center 3?>2 The 3D real space is divided into cubes whose size is 3 Å × 3 Å × 3 Å. 3?>2 is the center of a cube, which is introduced to specify its position. Note that "cube" and "zone" are different quantities: cube is defined in the 3D real space, and zone is done in the 3D RC space.
This quantity is cube-based bosentan-tail contact ratio, and shows spatial positions (cubes) where bosentan-N-terminal tail contact is formed frequently. When bosentan, whose centroid is involved in a cube centered at 3?>2 , contacts to the N-terminal tail by %, >> slice and ∆ F Axis of >> is divided into thirteen intervals (see Table S3). Volumes involved in the intervals are shown in Figure S13. Because the slices look like slices as in Figure S13, we call the interval "slice" or " >> slice" or "slice". The -th slice is denoted as ∆ F . ↓ and ← Vectors ↓ and ← are molecular orientation vectors introduced to define molecular orientation of bosentan in each snapshot. The vectors are defined on the core region of bosentan ( Figure 3a). Because the core region can be regarded as a rigid body approximately, the vectors are good measures to define the molecular orientation of bosentan. The two vectors are approximately perpendicular to each other, and the vectors are normalized as Thermal average of a quantity Ω in a cube centered at 3?>2 . This quantity shows spacial patterns of Ω in the 3D real space. Ω is defined arbitrarily by a user. See subsection 3.4 of SI.
Thermal average of ← and ↓ , respectively, in a cube centered at 3?>2 . This is an Tis quantity is denoted as residue-based bosentan-tail contact ratio, which expresses a ratio of bosentan contacting to residue of the N-terminal in the >> slice Δ F . See section 9 of SI.

39B K # ( , )
This quantity is denoted as atom-residue contact ratio between bosentan and N-terminal tail, which represents the ratio of a bosentan's heavy atom contacting to residue of the Nterminal tail in the slice Δ F . See section 9 of SI. See Figure S22a.

Section 1. Definition of a reaction coordinate (RC)
Consider two atom groups ! " and ! # (ℎ = , , ) in a molecular system. A single reaction coordinate (RC) (!) is defined by the distance between the centroids of ! " and ! # ( Figure S1). Superscripts and indicate simply that the two atom groups are pairing to define the single RC, and then, one can exchange the superscripts as: ! " → ! # and ! # → ! " without changing the value of (!) . The actual set of the atom groups for the present system is explained later.
Although we used centroids of atom groups to define RCs in this paper, RCs can be defined arbitrarily by the system's coordinates in theory. To define RCs, for instance, one can use angles that can be described by system's coordinates.

Figure S2
Figure S2: Generation of molecular system. A solid-line rectangle represents system at each stage, and solid-line arrow represents conversion of the system from a stage to the next stage. "Complex structure in figure S4", which is close to the X-ray structure of bosentan-hETB complex (PDB ID: 5xpr), may be called the "native-complex structure" in this paper if necessarily. Broken-line rectangles and broken-line arrows indicate actions applied to the system to proceed the stage. "Diffused conformations in figure S6" are initial conformations for GA-mD-VdMD simulations.

Section 2. Generation of the molecular system
To generate the bosentan-hETB molecular system, we can use either of two crystallographic complex structures 2 : Bosentan-hETB (PDB ID: 5xpr) and K8794-hETB (PDB ID: 5x93). The chemical compound K8794 is an analog of bosentan. Figure S3a illustrates the two complex structures similar to each other. Both compounds bind to the pocket of hETB, and importantly, their positions in the complexes are almost identical ( Figure S3b). K8794 has a long sidechain indicated by a magenta-colored broken-line circle in Figure S3b, which does not exist in bosentan. We select the K8794-hETB (PDB ID 5x93) complex to generate the computed system because 5x93 has a better resolution than bosentan-hETB (PDB ID 5xpr) complex: 2.2 Å and 3.6 Å for 5x93 and 5xpr, respectively. Then, we replaced K8794 of 5x93 by bosentan. The wild-type hETB consists of 442 amino-acid residues, which are numbered from 1 to 442 in UniPlot (https://www.uniprot.org/uniprot/P24530.fasta). In the crystallographic structure, the intracellular loop 3 (residues 304-310) of hETB was replaced by T4 lysozyme. Then, to generate the genuine hETB molecule, we removed T4 lysozyme and modeled the loop 3. In the crystallography, seven amino-acids of hETB were conducted to increase the thermostability of hETB (R124T, D154A, K270A, S342A, and I381A) or to prevent heterogenous palmitoylation in hETB (C396A and C400A). Then, we back-mutated those amino-acids to the wild-type ones. Besides, some missing atoms in the crystallographic structure were modeled. Although hETB is originally comprised of residues 1-442, the N-terminal residues 1-84 and the C-terminal residues 405-442 are missing in the X-ray structure of 5x93. In our previous simulation of the ET1-hETB system 1 , we treated residues 85-403 explicitly and found that the N-terminal tail (residues 85-89) of hETB captures ET1 prior to binding of ET1 to the pocket of hETB (the fly-casting mechanism 3,4,5 ). In the present study, we extended the N-terminal tail up to the 80-th residue to study the fly-casting mechanism further. Therefore, the N-terminal tail consists of residues 80-89 in the current study. We did not extend the C-terminal tail of hETB because the C-terminal tail is in the cytoplasmic side (no direct interaction to bosentan). Finally, hETB consists of residues 80-403 in the simulation. The N-and C-terminals of hETB were capped by the acetyl and N-methyl groups.
According to a method explained in our previous paper (Section 8 of SI of ref. 1), we embedded hETB in a membrane (the popc-lipid bilayer) and put four cholesterol molecules in the hETB-membrane interface because these cholesterols play an important role to stabilize the GPCR (i.e., hETB) tertiary structure. 6 Up to here, the system consists of bosentan, hETB, and membrane (cholesterol and POPC lipid molecules). The system generated above was immersed in a periodic-boundary box filled by water molecules. The sides of the box were parallel to the x-, y-, or z-coordinate axis (the x-y plain was parallel to the membrane surface), where the box size was 71.33 Å (x-axis), 71.33 Å (y-axis), and 132.17 Å (z-axis). Then, we introduced 28 sodium and 43 chlorine ions with replacing water molecules randomly by the ions to set the net charge of the entire system to zero and to set the buffer ionic concentration to a physiological one. The number of the final constituent atoms of the system is 69,062 (numbers of protein and bosentan atoms are 5289 and 68, respectively; number of cholesterols, POPC, and water molecules are 4 (296 atoms), 127 (17018 atoms), and 15440 (46320 atoms), respectively; number of sodium and chlorine ions are 28 and 43, respectively). This conformation is denoted as "Complex structure in periodic box" in Figure S2.
After energy minimization of the above system, a short NVT (constant volume and temperature) simulation was performed. Then, an NPT (constant pressure of 1 atm and temperature of 300 K) simulation was performed to equilibrate the box size. The resultant box size was 69.37Å × 69.37Å × 140.29Å. Figure S4 illustrates the resultant conformation, which shows that bosentan is close to the natively bound position because the generation of the molecular system started from the X-ray structure of 5x93. We denote this conformation "Complex structure of Figure S4 (native-complex structure)" in Figure S2. We emphasize that this complex is not the initial conformation of GA-mD-VcMD: As explained in the main text, two more stags are processed to obtain the initial conformation, where bosentan was completely separated from the binding pocket of hETB.
Last in this section, we note that the conformation of the N-terminal tail in the native complex ( Figure S4) should be regarded as an instantaneous conformation because the tail is an intrinsically disordered segment, whereas the other part of the hETB is determined in the crystal structure 2 . We report in this paper that the N-terminal tail fluctuates largely during the simulation. This result agrees with the fact that the Nterminal tail is missing in the crystal structure.  Figure S2, where this conformation is named "Complex structure of Figure S4 (native-complex structure)".

Subsection 3.1. Unified coordinate system and spatial density of bosentan around hETB
The present GA-mD-VcMD simulation produced 114 × 10 & snapshots as mentioned in the main text. To analyze the molecular structure and molecular motions from the snapshots, a unified coordinate system should be defined. Because the trans-membrane (TM) helix regions of hETB were structurally conserved well during the simulation, we superimposed the TM helix regions ( Figure S7a) of each snapshot on those of the native complex structure ( Figure S4). The other portions of the system (regions other than the TM helix regions of hETB, bosentan, cholesterols, POPCs, water molecules and ions) in the snapshots were moved according to this superposition. This procedure corresponds to define a unified coordinate system on the TM helix regions (a body-fixed coordinate system). This superposition was applied to all the snapshots, and physical quantities were calculated using those superposed snapshots. Benefit of GA-mD-VcMD is that an equilibrated thermodynamic weight (statistical weight) at a simulation temperature (300 K in the present study) is assigned to each snapshot (equation 31 of Ref. 7). Here, we denote the thermodynamic weight assigned to snapshot as ' . The position of the bosentan's centroid ( ) is referred to as ()* , and ()* of snapshot as ()*,' . A method to calculate a spatial density of ()* around hETB is given as follows: First, the 3D real space is divided into cubes whose volume ∆ is 3 Å × 3 Å × 3 Å. The position of a cube is represented by its cube-center ,-*. in 3D real space. Then, the snapshot is assigned to a cube that involves ()*,' . Assigning all the snapshots to cubes, we calculate the spatial density ()* ( ,-*. ) of the bosentan's centroid as: (S1) The function ()* ( ,-*. ; ) is a delta function defined as: Last, we normalized ()* ( ,-*. ) so that the maximum density is 1: Also, we define the spatial density of the centroid of a molecular portion other than the whole of bosentan. Denoting the centroid of of snapshot as ()2,' , the spatial density of is defined as: where ()2 is a delta function defined as: For instance, in this study, is set to the tip (Ser 80) of the N-terminal tail of hETB ( = ) or the tip (Lys 248) of the -hairpin of hETB ( = ℎ). Then, ()2 ( ,-*. ) is denoted as ()34 ( ,-*. ) and ()5! ( ,-*. ), respectively.

Subsection 3.2. Distance distribution function and radial distribution function
Here, we present a method to calculate a distance distribution function 667 ( ) and radial distribution function 867 ( ), where is a variable to define the functions. In this study, is the root-mean-square-deviation ] @/B , where ' and ' =.A are the -th heavy atom of bosentan in a snapshot and the native complex structure, respectively. For = , the summation is taken over all heavy atoms in bosentan and 04</ is the number of the heavy atoms. For = , the summation is taken over those in the bosemtan's core region (see Figure 3a in the main text for the definition of the core region), and 04</ is the number of heavy atoms in the core region.
The quantity 667 ( ) is a probability to detect the system in a distance window We note that the trans-membrane helices of the snapshot are superposed on those of the native complex structure in advance (see Subsection 3.1 of SI and Figure S7a) and that no further superposition is applied to the coordinates of the snapshot to calculate . Therefore, is contributed by three components: Translation of bosentan's centroid (or centroid of the core region) from the snapshot to the native complex, rotation around the bosentan's centroid between the snapshot and the native complex after the translation, and the conformational difference of bosentan between the snapshot and the native complex after the rotation. We note that converges to the translation with bosentan with bosentan going away from the native-complex position, and that the translation is regarded as a distance between two centroids between the snapshot and the native complex, which is denoted as ** in Figure S7c: 667 ( ) → 667 ( ** ) and 867 ( ) → 867 ( ** ).

Subsection 3.3. Bosentan-tail contact ratio
A quantity named "cube-based bosentan-tail contact ratio" is introduced as follows: We calculated the minimum heavy-atomic distance, *?3 , from bosentan to the N-terminal tail (Ace-Ser 80-Pro 81-Pro 82-Arg 83-Thr 84-Ile 85-Ser 86-Pro 87-Pro 88-Pro 89; Ace is the acetyl group introduced to cap the N-terminal of the N-terminal tail) for snapshot . If *?3 ≤ 5 Å, the substantial space between bosentan and the N-terminal tail is equal to or smaller than 1 Å (= 5 Å − 2 Å − 2 Å) assuming that the radius of a heavy atom is 2 Å approximately. This space of 1 Å is smaller than the diameter of a water molecule (3 Å). Thus, we judged that bosentan and the N-terminal tail were contacting substantially in the snapshot. Next, we assigned bosentan of the snapshot to a cube ,-*. that involved the centroid ()2,' of bosentan. Then, we defined the cube-based bosentantail contact ratio: where The ratio ̅ *?3 ( ,-*. ) quantifies the bosentan-tail contact ratio in the cube at ,-*. . The larger the ̅ *?3 ( ,-*. ) in the cube ,-*. , the higher the contact ratio in the cube. The maximum of ̅ *?3 is 1.0 (bosentan contacts always to the N-terminal tail when bosentan is in the cube), and the minimum is 0 (always no contact).

Subsection 3.4. Spatial patterns of quantity
Here, we defined a method to calculate spatial patterns of a quantity Ω defined by the coordinates of the system (i.e., Ω is a function expressed by the system's coordinates): Ω = Ω(x @ , y @ , z @ , … , x 3 dee , y 3 dee , z 3 dee ) , where x J , y J , and z J are, respectively, the x-, y-, and z-coordinates of atom in the system, and 099 is the number of constituent atoms in the system ( 099 = 69,062 in the present system). The quantity Ω is not necessarily defined by all the coordinates but may be done by some of the coordinates. The spatial density is defined as: where Ω ' is the value of Ω of snapshot . If the quantity Ω is a vector, then Ω is replaced by in the above equation.  Figure S8 demonstrates the resultant conformational distribution ,0L< v (>) , (5) , (M) w in the 3D-RC space, which was calculated straightforwardly from the simulations. The value of ,0L< is provided in a common logarithm. The highest density region (i.e., the lowest free-energy basin) is shown by the red-colored contours in the Figure: log @N [ ,0L< ] = −0.5, whose free energy is − ln[ ,0L< ] = 0.68 kcal / mol ( is the gas constant and = 300 K) measured from the lowest free-energy site. Note that the native-complex structure (small black sphere labeled L04':. in the Figure) is located at the periphery of the lowest free-energy basin. The blue-colored region (log @N [ ,0L< ] = −1.0; free energy of 1.36 kcal / mol) surrounds the lowest free-energy basin and it is stretched to the direction parallel to the (>) -axis. Remember that (>) controls the gate width of the hETB's binding pocket (Figure 1b). Figure S8 indicates that the free-energy slope is gentle to the gate motions. With increasing (5) and (M) , ,0L< decreased (i.e., free energy increased).  Table S2). Then, we convert the indices [  Table S2.
We performed 55 iterations as mentioned in the main text. Figure S9 illustrates the 0,, (0.25) region at iterations 1, 5, 10, 20, 40, and 55. This figure indicates that 0,, (0.25) increased rapidly in the first ten iterations. In iteration 10 ( Figure S9c), 0,, (0.25) involved both the native complex structure ( L04':. ) and the completely dissociated conformations ( O.P0=04. ). Therefore, we may quit the simulation at iteration 10 to obtain a rough and overall feature of ,0L< . In fact, the overall shape of 0,, (0.25) did not vary largely in iterations 10-55. However, we continued the simulation up to iteration 55 to save more conformations for analysis. Figure S9f

Section 6. Comparison of the current and previous distributions
We compared the conformational distribution ,0L< v (>) , (5) , (M) w with that obtained from our previous study of the ET1-hETB system 1 . The ligand-receptor approaching/departing was regulated by a single RC, B , in the previous study, whereas it was done by two RCs, (5) and (M) , in the present study. The gate opening/closing of the hETB's binding pocket was controlled by a single RC in both studies: @ in the present study and (>) in the previous study. Thus, the conformational distribution was presented two-dimensionally by a function ,0L< B6 ( @ , B ) for the ET1-hETB system, although it was expressed by ,0L< ( (>) , (5) , (M) ) for the bosentan-hETB system. For the comparison, we contracted ,0L< ( (>) , (5) , (M) ) to a 2D form ,0L< B6 ( (>) , (5) ) or ,0L< B6 ( (>) , (M) ) simply as: ,0L< B6 v (>) , (5) Then, the computed 2D distributions are normalized so that the largest value of the distribution, [ ,0L< B6 ] /01 , is set to 1: [ ,0L< B6 ] /01 = 1. Then, the 2D distribution can be converted to a potential of mean force as: where = or . Normalization is applied to ,0L< B6 so that the lowest B6 is 0: Figures S10a and S10b present B6 v (>) , (5) w and B6 v (>) , (M) w, respectively. In the both panels, we found a low free-energy path denoted by an arrow labeled " @ ", which was parallel to the (>) -axis. Therefore, the gate opening/closing motions follow the low free-energy path with gentle upward/downward slopes. Arrows labeled " B " and " R " represent motions along which bosentan approaches/departs hETB. We note that the free-energy slope along arrow B is slightly gentler than that along arrow R . This seems natural because bosentan can approach/depart the binding pocket of hETB with a less cost when the gate is in an open form. This was also observed in the landscape ,0L< ( @ , B ) in the previous study for the ET1-hETB system ( Figure 4A in Ref. 12). However, the scales of the free-energy landscapes in the two studies differed largely: The slope in the current bosentan-hETB landscape was about one-tenth of that in the previous ET1-hETB landscape.

Section 7. Division of axis into slices
We introduced a quantity ** in Figure S7c, which is the distance from the bosentan's centroid, ()* =.A , in the native-like complex to that, ()*,' , in snapshot . Here, we divide the 3D real space into thirteen slices (Δ J ; = 1, … ,13) along the ** axis with slice thickness of 5 Å. The actual slice ranges are listed in Table S3. Figure S13 demonstrates distribution of the bosentan's centroid ()* in each slice. Apparently, the distribution of the centroid in slices Δ @ ,…, Δ & are narrow, which indicates that those slices are involved in the binding pocket of hETB. In contrast, the distribution is wide in the slices of Δ S ,…, Δ @R . We note that the distribution of centroid varies suddenly at the boundary between Δ & and Δ S , and that this boundary corresponds well to the gate of the binding pocket, which was defined as the boundary between the blue-colored contour region ( ()* ( ,-*. ) ≥ 0.1) and the magenta-colored contour region ( ()* ( ,-*. ) ≥ 0.01) in Figure 2. Table S3.

Figure S13
Figure S13: Colored dot represent centroids of bosentan in snapshots. Color of a dot corresponds to a :: slice Δ 3 ( = 1, … ,13), to which the centroid belongs. See Table S3 for actual range of Δ 3 . For instance, Δ K indicated by broken line is in range of 20 Å ≤ :: < 25 Å. Black small sphere represents "centroid of bosentan in native complex" (cyan-colored stick model). Color assigned to Δ 3 is shown in inset, where four colors appear repeatedly as red → yellow → magenta → green with shifting slice as Δ 3 → Δ 3L9 . Five hundred centroids are shown in each slice, which are picked randomly from snapshots in the slice. The shown molecular structure is native complex.

Section 8. Conformations picked from some slices
Figure S14 exemplifies snapshots taken from Δ T (40 Å ≤ ** < 45 Å), and this slice is involved in the region of ̅ *?3 ≈ 0.1 of Figure 4. In Figures S14a and b, bosentan contacts to the tip of the N-terminal tail, while it does not in Figure S14c and d. As shown in Figure S14e, there are snapshots in that bosentan contacts to membrane without touching the N-terminal tail. We also found snapshots in that bosentan contacts to both the N-terminal tail and the membrane (Figure not shown). We do not show snapshots in Δ @N − Δ @R because the bosentan-tail contact was formed rarely.  Figure S15 illustrates snapshots in Δ U (30 Å ≤ ** < 35 Å), and this range involves the region of ̅ *?3 ≈ 0.7. In Figures S15a and b, bosentan contacted to the tail with more area than that in Δ T (Figures S14a and b). Therefore, the bosentan-tail contact becomes tightly with bosentan approaching the gate of the pocket. Figure S15c illustrates bosentan floating in solvent without contacting to the N-terminal tail, and Figure S15d does bosentan contacting to both the N-terminal tail and membrane. In Figure S15e, bosentan buried somewhat in the membrane without contacting to the N-terminal tail.  Figure S16 are those taken from Δ S (20 Å ≤ ** < 25 Å), and this slice involves the region of ̅ *?3 ≈ 0.9. This figure exemplifies that bosentan passes the gate with multiple ways. In Figures S16a and b, bosentan is situated in the white solidline rectangle of Figure S7b and sandwiched by the N-terminal tail and the -hairpin. On the other hand, in Figure S16c and d, bosentan is in the white broken-line rectangle of Figure S7b when passing the gate. Bosentan in Figure S16c contacts to the N-terminal tail, although bosentan in Figure S16d does not. Figure S16e is an example that bosentan cannot pass the gate because bosentan is blocked by the -hairpin and bosentan should take a roundabout route to reach the gate.   Figure  10A does not contact to bosentan, although the tail in Figures S17b and c contacts to bosentan. On the other hand, in Figures S17d and e, bosentan is passing the gate without contacts to either the N-terminal tail or the -hairpin. The bosentan-tail contact ratio decreases in the binding pocket as shown in Figure S14. This may suggest that that the role of the bosentan-tail contact decreases when bosentan is in the binding pocket.

Section 9. Residue-based and atom-residue bosentan-tail contact
To analyze the contact of bosentan to a residue of the N-terminal tail as a function of the ** slice Δ J , we introduce a "residue-based bosentan-tail contact ratio" ,L4 V g as follows: First, the minimum heavy-atomic distance between bosentan and residue of the tail is calculated for snapshot . Then, if the distance is smaller than or equal to 5 Å, then we judge that the -th residue contacts to bosentan. The residue-based contact ratio is defined as: where the summations in equation S11 are taken over all sampled snapshots. Remember that / is the thermodynamic weight assigned to snapshot . ( , ) is a delta function: The function / ( ) is defined as: Because the normalization of equation S11 was done by the entire probability (i.e., denominator ∑ / / of equation S11 is the probability contributed by all snapshots), we can compare the ratios between different slices: I.e., if ,L4 V g ( ) > ,L4 V gi ( ), the contact formed at residue in Δ J is stabler than that in Δ JX . Then, we convert this contact ratio to a potential of mean force ( ) for a freeenergy landscape as: Then, the inequality ,L4 V g ( ) > ,L4 V gi ( ) is converted to V g ( ) < V gi ( ). This quantity is again normalized so that the lowest value of is set to zero in a 2D free-energy landscape constructed by and Δ J : V g ( )Ž /'L = 0.
Next, to elucidate the bosentan's site specificity contacting to the N-terminal tail, we introduce an "atom-residue contact ratio" ,L4 V g ( , ), which is the contact ratio of a heavy atom of bosentan to residue of the N-terminal tail in the slice Δ J : where / ( , ) = F 1 (if residue contacts bosentan X s atom in snapshot ) 0 (else) . (S16) Because the normalization in equation S15 is done over snapshots involved in Δ J by the effect of (Δ J , ), the atom-residue contact ratio ,L4 V g ( , ) should not be compared among different slices. I.e., ,L4 V g ( , ) is designed to discuss the bosentan's site specificity in each slice.   Section 10. Inter-molecular native contacts in three X-ray structures. Figure S21a displays the X-ray structures of the bosentan-hETB and bosentan-hETB complexes. The K8794-hETB complex is not shown because it is very similar to the bosentan-hETB complex ( Figure S3b). Bosentan is located at a deeper position of the binding pocket of HETB than the main body (helical part) of ET1 is done ( Figure S21a). This may be an example that the ligands can bind to a deep or shallow position of the binding pocket of a GPCR depending on the ligand type 11 . However, all the three ligands adopt the same binding scheme as discussed in Ref. 2. Figure S21b clarifies that the sulphonyl group (SOO) of the sulfonamide of bosentan and K7894 overlap well with the carboxyl group (COO ? ) of the C-terminal residue (Trp 21) of ET1. Besides, the bulky sidechains (ring C defined in Figure S21c) of bosentan and K7894 overlap well with the sidechain of Trp 21. Therefore, the sulphonyl group and ring C of bosentan play important roles for stabilizing the complex structure.  Figure S22a illustrates the X-ray structure of bosentan-hETB focusing on the two oxygen atoms of the sulphonyl group of bosentan, which interact with the sidechain nitrogen atoms of three charged residues of hETB: Y of Lys 182 (3x33), Y of Lys 273 (5x39), and Z ′ of Arg 343 (6x55). The nitrogen atom ( 14) of the sulfonamide is also close to Y of Arg 182. We observed similar interaction patterns in the K7894-hETB complex (figure not shown). Figure S22b, which focuses on the ET1's C-terminal (Trp 21) of the ET1-hETB complex, demonstrates the interactions between the carboxyl group of Trp 21 and the sidechain nitrogen atoms of Lys 182, Lys 273, and Arg 343. Note that an atom corresponding to N14 of bosentan does not exist in Trp21 of ET1. N14 does not exist in the ET1's C-terminal. Distances among those atoms in the X-ray structures are shown in Table   S4.  ";$ is atomic distance between heavy atoms and , which belong to ligand (bosentan, K8794, or ET1) and hETB, respectively: = 14 or SOO of bosentan's sulfonamide as well as = Lys182 N j , Lys273 N j , or Arg343 N k . Exact definition for ";$ is presented in main text.

Section 11. Contacts between heavy atoms and
Given a snapshot , we define a delta function regarding the distance >;5 between atoms and : vSOO; Lys343N f w . The minimum and maximum of 〈 ,L4 ( ** )〉 are 0 and 4, respectively. The standard deviation of the number of contacts is calculated as:

Section 12. Formation of each native contact
To investigate the formation of each native contact, we introduced a function *O ( >;5 ; ( ** )), which is the contact probability as a function of >;5 using snapshots in a range ( ** ). Figure S23 presents six panels of *O ( >;5 ; ( ** )), and the panels are arranged so that bosentan is approaching the native-complex position from the gate of binding pocket: ** = 9.25 Å, 6.25 Å, 3.25 Å, 2.25 Å, 1.25 Å, and 0.25 Å. At ** = 9.25 Å, the probabilities for all the >;5 distances were considerably low ( Figure  S23a). Decreasing ** , small peaks appeared at >;5 ≈ 3 Å, which indicates formation of the atomic contacts: A peak for >;5 =`a a;\]^BUR[ o in v6.25 Åw (Figure S23b Figure S23d). In v1.25 Åw, the number of peaks increases although their peak heights were still low. Last in v1.25 Åw, all the four distances inhibited peaks at >;5 ≈ 3 Å or 4 Å, and the peak heights increased (Figure 16e). This suggests that the complex structure was stabilized by the inter-molecular interactions (i.e., enthalpy) at the bottom of the binding pocket as observed in the X-ray structure.